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Abstract 

We study the thermodynamic properties of an ideal gas of fermions in a harmonic oscillator 
confining potential. The analogy between this problem and the de Haas-van Alphen effect is 
discussed and used to obtain analytical results for the chemical potential and specific heat in 
the case of an isotropic potential. Step-like behaviour in the chemical potential, first noted in 
numerical studies, is obtained analytically and shown to result in an oscillatory behaviour of the 
specific heat when the particle number is varied. The origin of these oscillations is that part of the 
thermodynamic potential responsible for the de Haas- van Alphen effect. At low temperatures we 
show that there are significant deviations in the specific heat from the expected linear temperature 
dependence again as a consequence of the de Haas- van Alphen part of the thermodynamic potential. 
Results are given for one, two, and three spatial dimensions. 
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I. INTRODUCTION 

The experimental realization Q| of Bose-Einstein condensation in confined gases of atoms 
at low temperatures has been the stimulus for a wide range of experimental and theoretical 
investigations. (See Ref. jj] for a review.) By ingenious experimental techniques it is possible 
to prepare the atoms as either bosons or fermions. There have now been a number of 
experiments on trapped gases of fermions, as well as mixtures of bosons and fermions jjj. 

For a gas of fermions, it is a good approximation to model the system as a collection 
of non-interacting particles obeying Fermi-Dirac statistics confined by a simple harmonic 
oscillator potential. Unlike the case of bosons, the dominant s-wave scattering channel is 
suppressed making the effects of interactions less important 4j. With this simple theoretical 
model for the system, it is easy to obtain exact and simple expressions for the single- 
particle energy levels using the standard quantum mechanical result for the simple harmonic 
oscillator. This can serve as a starting point for a theoretical analysis of the thermodynamic 
behaviour of the system. 

Because the energy levels for a set of particles (bosons or fermions) confined by a simple 
harmonic oscillator potential form a discrete set, in order to calculate the thermodynamic 
properties we must be able to perform sums over the quantum numbers labelling the states; 
this can lead to difficulties in evaluation. By restricting the parameters of the system 
(temperature and oscillator frequencies for example) it may be possible to argue that the 
sums can be approximated with integrals, an approximation that renders the computation 
of analytical results easier. We will refer to the approximation of sums with integrals the 
continuum approximation, since it is equivalent to regarding the discrete energy spectrum 

as continuous. This approximation was first applied to the Fermi gas in the now classic 

pi 

paper of Butts and Rokhsar [4] . Even when one goes beyond this approximation, as we will 
show below, the continuum approximation can still be very accurate in some temperature 
regimes, and can capture the leading order behaviour if the temperature is not too low. 

The first study that showed there might be more features present bevond what is obtained 
by the continuum approximation was given by Schneider and Wallis [5|. These authors com- 
puted the thermodynamic expressions numerically, and their results showed that there were 
a number of step-like features present, not predicted by the continuum approximation. In 
particular, the chemical potential, which is a smooth function in the continuum approx- 
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imation, has a series of step-like jumps when a numerical analysis of the exact results is 
performed. This in turn can lead to step-like, or oscillatory, behaviour of other thermody- 
namic properties, such as the specific heat. Another important consequence of the work 
presented in Ref. j^J is that at sufficiently low temperatures there are significant deviations 
in the specific heat from the linear temperature dependence predicted jj| by the continuum 
approximation. The main purpose of the present paper is to analyze the situation analyti- 
cally to obtain the step-like behaviour of the chemical potential, and to apply this method 
to the evaluation of the specific heat. 

In add.t.n t o the two stud.es ahead, onioned QQ, there have been a number of other 

n n 

related theoretical works. Ref. [6| generalized the calculations of Ref. [5|| to a gas with two 
spin states with an interaction between the two states. (We will make some brief comments 
on the effect, of interactions in Sec. ED A number of authors Q fl fl fl Q Q Q 
have obtained analytical approximations in a number of cases using path integral, Green's 
function, or density matrix methods. 

In the present paper we will concentrate on the behaviour of the chemical potential and 
specific heat, and will use the evaluation of the thermodynamic potential as the basic starting 
point. The key observation is that for a simple harmonic oscillator confining potential, the 
evaluation of the thermodynamic potential is mathematically very similar to that of a gas 



of charged fermions in a constant magnetic field. Landau [1J| was the first to show that 
the thermodynamic properties of electrons in a constant magnetic field undergo oscillations 
whose period is determined by the inverse of the magnetic field strength. This gave a 
theoretical basis for the observations of de Haas and van Alphen, now referred to as the 
de Haas- van Alphen effect. (See Ref. for a good discussion of the de Haas- van Alphen 
effect and the important role it plays in condensed matter physics.) What we will show 
here is that for the trapped Fermi gas the step-like behaviour found in Ref. |5] is completely 
analogous to the de Haas- van Alphen effect. Just as in the de Haas- van Alphen effect, where 
use of the continuum approximation would miss the oscillations, here too we must go beyond 
the continuum approximation. 

The outline of our paper is as follows. In Sec.|H]we describe the general method following 
the classic analysis of Sondheimer and Wilson lira. There are some differences with the 

n 

de Haas-van Alphen case considered in Ref. 16] that we will describe. We will show how 
the continuum approximation comes about in the method, and even find the next order 
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correction to it, beyond the leading order result given in Ref. jj|. In Sec. II I II we analyze 
the de Haas-van Alphen contribution to the thermodynamic potential for a 3-dimensional 
isotropic harmonic oscillator potential, and obtain approximate analytical results for the 
chemical potential. This is then used to examine the specific heat. Results for the 1- and 
2-dimensional gases are given in Sec. I1VI and the results compared with the 3-dimensional 
case. Finally, Sec. IS contains a brief summary and a short discussion of the main results. 



II. GENERAL METHOD 



We will begin with the grand canonical ensemble and the thermodynamic potential f2 
defined by 

n = -Tj2f(E n ), (2.1) 

n 

where 

/(£)=ln[l + exp(-/3(£-/i))] . (2.2) 

We use the usual notation {3 = 1/T, with T the temperature, and choose units such that 
the Boltzmann constant is equal to one. Take the Laplace transform of f{E), and call it 
(p({3), so that 

oo 

<p(/3) = I dEe^ E f(E). (2.3) 
o 

The inverse Laplace transform of Eq. ()2.3|) reads 

c+ioo 

f(E) = ^- J dp<J* B <p{p), (2.4) 

c— ioo 

where c is an arbitrary constant chosen so that the integration path lies to the right of any 
singularities of ip(/3). We can now use Eq. ()2.4|) in Eq. ()2.1j) . The sum over the energy levels 
that results can be related to the partition function for the /i = system, defined by 

Z(f3)=Y, e ~ PEn ■ ( 2 - 5 ) 

n 

We can then write fl defined by Eq. (j2.1J) as 

c+ioo 

Q = ~ f d(3Z(-(3M(3) . (2.6) 

c—ioo 



It is necessary here to regard Z(/3) as a function of j3 with j3 viewed as a complex variable. 
Z(—(3) is the result obtained from this complex function by analytic continuation of the 
definition Eq. (J23D to Sfc(/3) < 0. 

We now define 3(E) to be the Laplace transform of f3~ 2 Z((3). (The factor of f3~ 2 is for 
later convenience |l6|.) This results in the definition 

oo 

Z((3) = {3 2 J dEe-? E 3(E). (2.7) 
o 

The inverse Laplace transform gives 

c+ioo 

3(E) = hi I d P ePE p~ 2z w- ( 2 - 8 ) 

c— ioo 

Substitution of Eq. (|2.7jl into Eq. ()2.6j) followed by the use of Eq. (12. 4j) results in 

oo 

Q=-T [ dEl{E)^0-. (2.9) 



o 



The Fermi-Dirac distribution function F(E) is 



F(E) = + (2.10) 

dE 



-T^-f(E) , (2.11) 



if we use Eq. (|2.2j) . This leads to 



fi = / dE3(E)^-F(E) , (2.12) 



o 



<9£ 

giving the key starting point in the Sondheimer- Wilson 16] analysis of the de Haas- 
van Alphen effect. The method rests on an evaluation of 3(E) defined by Eq. ()2.7|) or 
Eq. (j211) and its use in Eq. (J2~T21) . 

To evaluate 3(E) in Eq. (J2.8)) we require the partition function, or at least some infor- 
mation about the singularities of (3~ 2 Z(f3) in the complex /3-plane. For the case of a simple 
harmonic oscillator potential this information is very easy to obtain. We will consider the 
D-dimensional oscillator potential 

1 D 

l/(x) = -m J>fx 2 , ( 2 - 13 ) 

3=1 



with x = (xi, . . . ,xd) the D spatial coordinates. The energy levels are simply given by (in 
h — 1 units) 



D 

E- 



E(% + ^W (2-14) 

3=1 V J 



where n = (ni, . . . , no) and rij = 0, 1, 2, . . . for j = 1, . . . , D. We then have Z(f3) in Eq. (|2.5|) 
expressed as a product of geometric series that are easily summed to give 

^)=ri(T3^)- <"5) 

3=1 ' 

Because of the neglect of interactions, and the simple form of the potential, it has been 
possible to evaluate the partition function in closed form. 

To obtain 3{E) we need to know the singularity structure of j3~ 2 Z(/3). It is now obvious 
from Eq. (|2.15|) that f3~ 2 Z{(3) is a meromorphic function. There is a pole at (3 = of 
order D + 2 as well as a series of poles along the positive and negative imaginary /3-axis at 
P — ftkj = 2nikj/ujj with kj = ±1, ±2, . . .. The order of these poles depends on the relative 
ratios of the harmonic oscillator frequencies. If none of the oscillator frequencies are rational 
multiples of any of the others, all of the poles, apart from the one at (3 — 0, will be simple. 
However, if some of the oscillator frequencies are rational multiples of the others, some of the 
poles away from the origin will be higher order. This makes a determination of the residues 
of poles along the positive and negative imaginary /5-axis for a general harmonic oscillator 
potential an inelegant treatment of special cases. We will return to this problem later. 

We can write 3(E) in Eq. ()2.8|) as 

3(£) = 3o(£) + 3r(£), (2.16) 

where we close the contour in the left hand side of the complex plane and use 3o 

(E) to 

denote the contribution to the integral from the pole at ft — 0, and 3r{E) the contribution 
coming from the rest of the poles along the positive and negative imaginary /5-axis. When 
Eq. (|2.16|) is used in Eq. (|2.12|) we obtain the thermodynamic potential as 

Q = Q + n r (2.17) 

in an obvious way. As in the de Haas- van Alphen effect, the oscillations in thermodynamic 
quantities will come from Q r , but we will first examine Qq. Q r will be evaluated in the next 
section. 
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To obtain 3o(E) we need the residue of [3~ 2 e^ E Z{f3) at (3 = 0. For general values of 
D and arbitrary frequencies u>j this is difficult to write down in any simple way. We will 
concentrate on just the cases D = 1, 2, 3 here, although there is no intrinsic reason why the 
method cannot be extended to other values. A straightforward calculation yields 

3o(£) = g->rD = l, (2.18) 
ME) = (4£ 2 - oof - uf) for D = 2, (2.19) 

3o(£) = 7—^ [240£ 4 - 120£ 2 H 2 + 4 + ul) + 7(cu 4 + u\ + cu 3 4 ) 

57bUa;ia;2a;3 L 

+10(w?w| + cj 2 cj 2 + w 2 ^) for D = 3. (2.20) 
These results can now be used in Eq. ()2.12|) to obtain 

oo 

n = y dE3 (E)^F(E). (2.21) 
o 

As T — > 0, meaning that — > oo, the Fermi-Dirac distribution function approaches a 



step-function, — > 6>(/i — -E), so that -^F(E) — > —5(E — jj). This crude approximation 



results in 

-3o(/i) , (2.22) 
and is really the first term in a systematic expansion due to Sommerfeld. (See Refs. 
Q] for three derivations of the Sommerfeld expansion.) In our case we find 

* -3o(A0 - y^V) - ^ 4 3£'V) . (2.23) 

Although the general Sommerfeld expansion contains higher order terms, these vanish here 
because "3o{E) is no more than quartic in E due to our restriction that D < 3. For arbitrary 
D, the expansion Eq. ()2.23|) will contain more terms with increasing powers of T. For low 
values of T these should be less important than those indicated in Eq. ()2.23|) in any case. 
Using Eqs. (j2. 18112. 20 jl in Eq. (ZZfl results in 



Q ~ -i-Gu 2 + yT 2 ) + |L , for JJ = 1 , (2.24) 
fi ~ + ^-^(^ 2 + 4 - 4vr 2 T 2 ) , for D = 2 , (2.25) 

bUiLU 2 2ALU1UJ2 

Q ~ ^ + ^ (a; 2 + ^ 2 + uj\ - 4tt 2 T 2 ) \l(uj\ + u4 + 

+ 10(cj 2 cj 2 + cj 2 cj 2 + ulul) - 40tt 2 T 2 (cj 2 + cj 2 + cj 2 ) + 112tt 4 T 4 ] , for D = 3 .(2.26) 
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The thermodynamic properties of the system follow from a knowledge of the thermody- 
namic potential. In particular, the average particle number N is given by 



(2.27) 



and the internal energy U is given by 

(2.28) 



In the most physically interesting case we take D = 3. (The cases of D = 1,2 will be given 
later in Sec. II VI ) If we temporarily ignore the contribution Q r to Q, we have N ~ N where 

*--(S)L> (2 - 29) 

- 7-^ --^ H 2 + uj\ + ujI - 4tt 2 T 2 ) . (2.30) 

OLO1UJ2UJ3 24:LJiUJ2td3 

For large values of fi, if we keep only the term in Eq. (|2.3(J|) of order /i 3 , we find (with 

/i 3 ~ Quo^uo^N , (2.31) 

showing that /i oc A^ 1 / 3 . For large values of N we will have /i much larger than the average 
oscillator frequency, consistent with the assumptions made in deriving these results. It is 
easy to obtain the next order correction to Eq. ()2.31j) . 

t r + M±d±Mzi^!>), (2 . 32) 



U1OJ2UJ3 \ 4(6A^a;iu; 2 u;3) 2 / 3 

showing that the correction to Eq. (|2.31|) becomes increasingly unimportant for large values 
of N. If we use Eq. ()2.32)1 and rewrite the result in terms of the Fermi energy defined by 

Ep = /i(T = 0) , (2.33) 

it is easy to see that 

/ 7T 2 T 2 \ 

This agrees with the result in Ref. j4| who used the continuum approximation and only 
the T-dependent part of the second term in Eq. (|2.30p . We have demonstrated here that 
Eq. (|2.34j) holds true even when the next order approximation for the density of states 
is included in the continuum approximation. However, as we will see in the next section 
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there are important corrections to this from terms that come from beyond the use of the 
continuum approximation and support the numerical investigations of Ref . \^ . 

We can also calculate the contribution to the internal energy, that we call Uq, coming 
from Q . Using Eq. (J2.26j) in Eq. (J2.28)) and then eliminating fi in favour of N with Eq. (J2.32j) 
results in the following approximation : 

U * h**w)V>N* + + A' 2 ' 3 • (2.35) 

assuming that N is large. (The intermediate result for Uo in terms of \x can be found in 
Eq. (|3.25|) below for the isotropic potential.) 
The specific heat can be found using 



(2.36) 



\dT J 

The contribution coming from the approximate result in Eq. (|2.35|) is then seen to be 

C ~ 7r 2 (6cui^ 3 )- 1/3 Ar 2 / 3 T . (2.37) 



As before this agrees with the result quoted in Ref. [4|. Again we will find that at low 
temperatures the part of Q that has not yet been considered, namely Q r , alters this expected 
behaviour, confirming the numerical results of Ref. Q|. We consider the evaluation of Q r 
and its effect in the next section. 



III. ISOTROPIC 3-DIMENSIONAL POTENTIAL 



A. Thermodynamic potential 

We will examine the case of the isotropic harmonic oscillator potential with 

^1=^2 = ^3 = ^, (3.1) 

in which case the results of the previous section simplify. Our aim here is to calculate fl r 
that has been neglected up to now from the analysis. To do this we need to evaluate the 
contribution to Q coming from the poles of j3~ 2 e^ E Z{(3) along the imaginary /3-axis away 
from (3 = 0. From Eq. ()2.15|) . because 



there will be poles of order 3 at (3 = (3k defined by 

(3 k = —k (3.3) 

for k = ±1, ±2, .... It is a straightforward matter to evaluate the residues of /3~ 2 e /3E Z((3) 
at (3 = f3k and show that their contribution to Eq. ()2.8j) gives 

ME) = £ (6 4- «V - cafrkE/v) + £ &(SriSM • 

fc=i ^ ' fe=i 

(3.4) 

It now remains to use this in Eq. (j2.12j) and to try to extract something useful from the 
result. 



We have 



(3u f^(-l) fe 4 (3 



64vr 4 ^ A; 4 8vr 3 ^ A; 3 

fc=l k=l 



where 



oo 



, , ,, ,,, 4?r 2 A: 2 ^ 2 \ cos(2nkE/uj) , . 

A k = dE[6 + n 2 k 2 — — ano/rn m » ( 3 - 6 ) 





oo 



cosh 2 [±/?(£-//)] 



5 fc = IdEE T 2 \T:J \t (3-7) 



sin(27r£; i?/u;) 

cosh 2 [±5(£-/i)] 
u 

In these last two integrals we can make the change of variable 

E = H + jO, (3-8) 

to new integration variable 6. The lower limits on the integrals in Eqs. ()3.6|) and ()3.7|) then 
become — /3/i/2. If we look at low enough temperatures, specifically T << /i as we have 
already assumed, then to a good approximation (up to exponentially small terms) we can 
replace the lower limits on the integrals defining A k and B k with — oo. The approximate 
results then become 



oo 

2 f dO 



(3 J cosh' 9 



2l2 4tt 2 A;V \Qix 2 k 2 ii6 lQii 2 k 2 e 2 
6 + ir 2 k 2 - 



lu 2 (3lu 2 (3 2 lu 2 

xcos - + -^) , 3.9 

2 f d9 ( 20 \ [2nku 47ik9\ 



OO 



f3 J cosh 2 6 V p J \ u f3uo J 
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The integrals in Eqs. ()3.9|) and ()3.10|) may be evaluated exactly using residues. All of the 
expressions required may be related to the basic integral 



cos(a6 + b) 
cosh 2 8 



d9 



ixa cos b 
sinh(|a) 



(3.11) 



and derivatives with respect to the parameters a and b. After some straightforward calcu- 
lation, using the results for Ak and Bk found from Eqs. ()3.9|) and ()3.1()|) as described, it can 
be shown that 



oo 



E 



-r 



Iti 2 (3 3 uj 2 ^-f k 3 sinh , , 

fc=l v (3u> 



8vr 4 A;^csch z 



27T 2 fc ' 

+2(3 2 u 2 + 7r 2 k 2 (A7r 2 + (3 2 u 2 - 4/3V) COs(27Tkfi/i0) 



2n 2 k 



47r 2 fc/3a;coth 



2n 2 k 



+AnkPjj (3u + 2tt k coth 



sin(27rkjj,/uj) 



(3.12) 



The presence of the trigonometric functions in Eq. ()3.12j) is responsible for the oscilla- 
tions that occur in thermodynamic quantities if we go beyond the leading order continuum 
approximation. The presence of the hyperbolic functions in Eq. ()3.12|) with argument 
means that unless f3u > 1 these oscillations will be suppressed. The oscillatory behaviour 
should show up for T < uj and become more prominent as T is reduced. 

The inclusion of Q r in the expression used for the thermodynamic potential will lead to 
corrections to the thermodynamic behaviour of the system beyond what is found using the 
continuum approximation of Sec. |H] We will look first at how the chemical potential is 
affected. 

B. Chemical potential 



We can define (from Eq. (j2~27j) ) 



djj, 



(3.13) 



and then use Eq. (|3.12|) to obtain 



N r 



Tt 



E- 



4/33(^ S inh (^) 

k=l v flu ' 



4/?V - 4tt 2 - (3 2 lu 2 - 87r 2 csch 2 



2n 2 k, 



sin(27r kfi/u) 
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27T 2 fc "I 

+87T/3/X coth (— — ) cos(2irk fi/u) \ . (3.14) 

This will make an oscillatory correction to the contribution Nq for the average number of 
particles found using the continuum approximation in Sec. ITT1 
The chemical potential may be found by solving 

N = N + N r (3.15) 

for fi, where No is given by (see Eq. ()2.3U|) with uji = uj 2 = oj^ = to) 

( " 6) 

Obviously the complicated dependence on /1 in N r as given in Eq. ()3.14|) renders an analytical 
evaluation of /x difficult. We can simplify by using the assumption that /1 is large ( since we 
have already assumed that /3/j 3> 1 and // ^> w), and keep only the leading term in /i. This 
gives 

N -(^ft smh(^) • (3 - 17) 

k=l \ pu> 1 

The continuum approximation in Eq. ()3.16|) still gives the leading contribution to iV 
for large /1; however, N r in Eq. ()3.17j) can be more important than the sub-leading term 
(proportional to fi) in Eq. (|3.1fij) if the temperature is low enough. For /3u of the order of 1 
or less, the sum in Eq. (|3.17j) is well approximated by simply the first term and it is clear 
that the oscillations, although present, will have a very small amplitude. We would therefore 
expect that for the temperature range T > u>, Eq. ()2.31|) or Eq. ()2.32|) would provide a good 
approximation for the chemical potential. However as the temperature is reduced, so that 
(3uj 3> 1, the amplitude of the oscillations coming from Eq. ()3.17j) become more pronounced 
and must be taken into effect. 

It is possible to find an asymptotic expansion for N r in Eq. ()3.17j) valid for j3u ^> 1 (or 
T <C uS). After some calculation, it can be shown that 

N r ^ £^ jtanh (^(2/2 - 1)) + 1 - 2/2 j , (3.18) 

where 

H= , (3.19) 
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with [x] denoting the largest integer whose value is less than or equal to x. (Thus < p, < 1 
can be assumed in Eq. (|3.18jl .) Similar asymptotic expansions can be found for the sub- 
leading terms in N r given in Eq. (|3.14j) ; however we will not give them here. It can be 
shown by numerically evaluating the sum in Eq. (J3.17|) and comparison with the analytical 
approximation in Eq. (|3.18|) that Eq. (J3.18|) does give an accurate result for large values of 

If we take the limit flu — > oo in Eq. (j3.18|) . the result simplifies further to give 

' 2 /V » , 1" 



2u 2 



cu + 2 



(3.20) 



Strictly speaking, this further approximation is only valid if + |] is not equal to an 
integer. In this special case, N r in Eq. (J3.17)) or Eq. (j3.18j) can be seen to vanish, and we 
must look at the sub- leading contributions to N r that follow from Eq. (J3.14)) . 

If we use Eq. (jgHH for N r in Eq. <$JB) along with Eq. (jHUBD for jV , it is possible to 
show that the solution for /i is given by the simple expression 

Z 1 r/cAr\l/3l , 1 



[(6Ar)V3] + (3 21) 

The presence of the greatest integer function in Eq. 1)3.21)1 leads to the step-like behaviour 

n 

first found in the numerical studies of Our results provide a confirmation of this be- 
haviour by analytical means. For large values of N, our result in Eq. ()3.21)) shows that these 
steps will occur roughly for N ~ £ 3 /6 where £ is an integer. This agrees with the "magic 
numbers" found by |5[ that occurred for N = £[£ + 1)(£ + 2)/6 if we make N, and hence £, 
large enough. We have therefore seen how the step-like behaviour of the chemical potential 
comes about in an analytical way, and traced its origin back to the same type of terms that 
are responsible for the de Haas- van Alphen effect. 

In Fig. ^we show the result for fi/oo plotted as a function of iV over a range of N. We 
have taken iV large, but the same type of behaviour can be found for smaller values as well. 
The continuum approximation for \i is shown as the smooth solid, almost straight, line, and 
the simple analytical approximation in Eq. 1)3.21)) is superimposed on it. The steps occur 
at the magic numbers as predicted. We have displaced the curves for different values of (5ui 
for clarity to show the trend towards the step-like behaviour as the temperature is reduced. 
As the temperature is increased, the amplitudes of the oscillations decreases as mentioned 
above. 
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FIG. 1: (color online) This plot shows (j,/lo plotted over a range of N for three sample temperatures. 
The solid, almost straight, line shows the result found using the continuum approximation for the 



found from the very low temperature approximation of Eq. 1)3.21(1 . The smooth sinuous curves give 
the results for T = 0.1u,T = 0.05a; and T = 0.01a; . The last two curves have been displaced from 
the first one for clarity. 

C. Specific heat 

We will first calculate that part of the energy that arises from Q r . Define 



particle number resulting in fi/cj 



(6./V) 1 / 3 . The step- function superimposed shows the result 



(3.22) 



and use Eq. (|3.12|) for Q r . After some calculation it can be shown that 
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2, ,2 



3tt^ 



(3 2 u 3 



E 



-l) fc cosh 



sinh 



2 /27r 2 fc^ 



COS 



/ 27r/c/A 



\ UJ 



(3.23) 



fc=l V /3aj / 

where we only include the two leading order terms in fj, since it is these two terms we will 
need to calculate the leading contribution to the specific heat. 

The specific heat was defined in Eq. (|2.36j) . It is straightforward to show that this 
expression is equivalent to 



(9N)\ 
\dTj\ w 







(SO 


7> 



(3.24) 



which is more useful for explicitly calculating C. The presence of the second term on the 
right hand side complicates the evaluation of C, but we will obtain an expansion for C in 
powers of \x. 

If we look at the first term on the right hand side of Eq. ()3.24|) and use U — U + U r with 
U r given by Eq. lET23|) and 

d 



d(3 

8cj 3 16cj 3 



An 2 T 2 ) + 



1 



we obtain 



1920cj 3 
vr 2 /i 2 T fdU, 



(17a; 4 + 40tt 2 uj 2 T 2 - 1 12vr 4 T 4 ) (3.25) 



fl,UJ 



2u 3 



+ 



dT 



(3.26) 



From Eq. ()3.23|) . counting powers of fi, it can be observed that the second term on the right 
hand side of Eq. ()3.26j) will contain expressions in /i 3 , fi 2 , . . .; thus, it appears as if the leading 
order behaviour of the specific heat will be /i 3 , rather than fi 2 as predicted by the continuum 
approximation. However, we will show that the /x 3 part of Eq. ()3.26|) cancels with a similar 
expression coming from the second term in Eq. ()3.24|) leaving the overall leading behaviour 
of the specific heat as /i 2 rather than /j, 3 . In any case, we must work consistently to order /i 2 , 
and to this order only the two contributions to U r that we have written down in Eq. ()3.23j) 
are necessary. We will drop all terms that result in explicit factors of fi in C that are of 
order /i and lower. 

Using Eq. (|3.25j) we have 



dU 
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+ 



dU r 
dfx 



(3.27) 



where the second term on the right hand side contains terms involving /i 3 ,/i 2 
Eq. (EHo) and Eq. (ETToT) we have 

ttVT fdN r 



fdN\ 
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fl,LU 



3u 3 
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From 



(3.28) 
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with the second term on the right hand side involving /i 2 ,/i, 

■ 2 'ONr 



as well as 



f dN 



T,t 



2u 3 



+ 



dfx 



(3.29) 



with the second term on the right hand side involving fi 2 , /i, . . .. A simple counting of powers 
of /j,, using the expressions for U r and N r , shows that the second term on the right hand side 
of Eq. ()3.24|) does begin at order /i 3 . 

It now remains to use the explicit results for U r given in Eq. ()3.23|) and N r given in 
Eq. (j3.14j) and evaluate C for large /i keeping terms of order /i 3 and /i 2 . After some calculation 
it can be shown that the order fj, 3 terms cancel leaving 



(3.30) 



with 



E, i cosh^i, ir 2 kT ( 1 2\1 / , £t\ 
-l fc <^ ^ — — + ^— \ cos (27rfc-) , 3.31 



where 
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OG 



E, i 1 2ix 2 kT cosYiOu 1 / , u\ , 
-1 — ^ Uin (2vrA;- ) , 3.32 
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k 

^ /3c t ;sinh6'fc 
k=i y K 



k=l 

00 - -\kA„2 



kb^r cos ( 2rt 5) • (3 ' 33) 



4 - ^ (3.34) 

has been defined to save a bit of writing. The continuum approximation of Sec. |H] can be 
regained by dropping all of the terms in £1,2.3 that have arisen from the de Haas-van Alphen 
part of the thermodynamic potential. This gives the familiar linear dependence on temper- 
ature P|. 

The numerical calculations of [5] showed that as the temperature got sufficiently small 
there was a significant departure from the linear temperature dependence in the specific 
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heat. Once again we will show that this follows from our analytical method, and that the 
origin of this behaviour is in the de Haas-van Alphen part of the thermodynamic potential. 
To do this we will evaluate the asymptotic expansion of the three sums defined in Eqs. (|3.31l - 
13.33)) for large values of f3uj. Leaving out the technical details of this for brevity, we find the 
approximate forms 



3(3 3 oo 3 (2fi - l) 2 
167r 2 cosh 2 (£* (2/2-1)) 



Si * — - 2 \T " x " 1 > ( 3 - 35 ) 



E 2 c , , (3.36) 

167rcosh 2 (^(2/2-1)) 1 ; 



go; 

4 cosh 2 (^(2/2-1)) 



S 3 ~ — Jp- , (3.37) 



with /2 defined as in Eq. (|3.l9)l . The results in Eqs. (j3.35fJ3~3"Tj) can be checked against a 
numerical evaluation of the sums defined by Eqs. (j3.3lH 3.32) and found to be accurate for 
(3uj ~ 10 and /2 not too close to or l. Once (3uj ~ 100 the results become very accurate 
even for JJl close to and l. Thus for T < u/100, the simple expressions in Eqs. (j3.35fJ3~3"Tj) 
become reliable approximations for £1,2,3- 

If we use Eqs. fl3.35fJ3~37]) in the expression for the specific heat in Eq. ()3.30|) the result can 
be shown to vanish. The de Haas-van Alphen contribution to the specific heat cancels the 
continuum approximation to the leading order we are working to. The specific heat therefore 
vanishes as T — > faster than T. This is completely consistent with the numerical results 
found in Because the de Haas- van Alphen approximation, as well as the asymptotic 
evaluation of the sums leading to Eqs. (|3.35Hl337j) neglect terms that are exponentially 
suppressed, we suspect that the specific heat vanishes like e~ au ^ T for some constant a as 
T — >• 0, but we have not been able to establish this in any simple way. Further support 
for this belief follows from the result we are able to establish for the 1-dimensional gas in 
Sec. IIVAI below. A more refined estimate of the sums, as well as the de Haas- van Alphen 
contribution would reveal the exact nature of the T — > limit. 

For T > cj/100, but still small, the results in Eqs. ()3.35Hl3"3T|) start to become less 
reliable. As a check on our results against the numerical ones of p| we plot the specific heat 
as found from Eq. ()3.3()|) to demonstrate the de Haas- van Alphen oscillations. This is shown 
in Fig. 121 As the temperature increases we do find, as expected, that the contribution from 
Ex 2 3 becomes smaller exponentially, and the linear behaviour with temperature is regained. 
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FIG. 2: (color online) This plot shows C/(6iV) 2/3 plotted over a range of N for four sample 
temperatures T = 0.1a;, T = 0.2a;, T = 0.3a; and T = 0.5a;. As the temperature increases the 
oscillation amplitude decreases and the curves approach the continuum limit of 1 in the scaled 
specific heat. As the temperature decreases there are significant deviations from the result found 
from using the continuum limit, and for very small temperatures the specific heat starts to become 
vanishingly small. 

IV. ONE AND TWO DIMENSIONS 

We will examine the thermodynamics of trapped Fermi gases in one and two spatial 
dimensions, since these cases may be of relevance as limiting cases of the 3-dimensional gas 
in some situations. Because the methods used are similar to those described above in the 
3-dimensional case we will be brief here and only exhibit key results that show a difference 
with results already obtained. 
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A. One dimension 



With D = 1, Z(f3) in Eq. (|2.15j) has simple poles at (3 = /3& with j3k still defined by 
Eq. ()3.3j) . It is easy to show that the contribution from the poles with k ^ to ~b{E) is 

fe=i v 7 

The contribution from the pole at /? = was given in Eq. (|2.18|) . After some calculation, 
making the approximation /i ^> T, it can be shown that 



1 (-l) fc cos(27rfc/i/q;) 

Pf-! fcsinh(^) ' l ' J 

k=l \ pw ' 



The calculation is very similar to the D = 3 case, so details will not be given here. The 
continuum approximation for Q, called Qq, was given in Eq. (|2.24|) . 

Using Eq. 1)2.24)1 and Eq. ()4.2jl in the general expression for the average particle number 
N in Eq. flZZZD results in 

at ~ ^ j_ 277 V" (-l) fc sin(27rfc/i/^) 

fe=i v ' 

The first term on the right hand side comes from Qq, and the second term from Q r , so that 
the continuum approximation would yield simply 

fi~uN (4.4) 

in this case. Again for large particle numbers we expect /i > w to be a valid approximation. 
For T <C w we can again find the asymptotic form of the sum in Eq. ()4.3)) and obtain 



N 



1 1 , {(3u (2u n 
- + - tanh I — M- - 2 



(4.5) 



2 2 I 4 V ^ 
and this proves to be very accurate for low values of T. 
The internal energy can be found to be 

with 

2tc ^ k (a sin(27rA;/i/a;) n cosh# fc . .1 
C/r ~ -r > (-1) < - — . : V + 7: o — sm(2vrA;u/a;) ^ (4.7) 
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and 6k denned by Eq. (|3.34|) . We can calculate the specific heat using Eq. ()3.24j) with the 
result 

The sums £1,2,3 are the same as those given in Eqs. (|3.31Hl3~34"j) . The result in Eq. (J4.8|) 
is very similar to that found in Eq. ()3.30|) in the 3-dimensional case apart from the overall 
factor that can be recognized as the continuum approximation for the specific heat of the 
1-dimensional gas. We can therefore conclude immediately that as T — > 0, the specific 
heat vanishes faster than the linear temperature dependence deduced from the continuum 
approximation. 




T/oj 



FIG. 3: (color online) This plot shows the specific heat plotted over a range of temperature for 
the 1-dimensional gas. As the temperature increases the specific heat approach the continuum 
limit that behaves linearly on the temperature. As the temperature decreases the result deviates 
substantially from the linear result and the curve decays exponentially fast as found in Eq. (|4.9|) . 

Because of the simplicity of the D = 1 case, we can shed some light on the behaviour of 
the specific heat as T — > 0. We note that the continuum approximation for the specific heat 
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in Eq. (|4.4|) is also a solution when the de Haas- van Alphen contributions to N are included, 
since the sin term vanishes for integral N. The sum £2 defined by Eq. ()3.32|) vanishes as 
well if /i = uN, and we can use the approximation for £1 given by Eq. (|3.35|) . valid for 
T <C uj to obtain from Eq. (J4.8|) the simple result 

This demonstrates that the specific heat does vanish exponentially fast as T — > 0, not 
linearly. The role of the de Haas- van Alphen part of the thermodynamic potential is again 
responsible for this behaviour. As a check on this conclusion we have plotted the specific 
heat as a function of temperature in Fig. |3] The results can be seen to be consistent with our 
analytical result in Eq. (|4.9|) as T — > 0. For larger values of the temperature the specific heat 
approaches the linear temperature dependence predicted by the continuum approximation. 
This can also be seen from the expression obtained in Eq. (|4.8jl since as f3us becomes small, 
the sums £1.2,3 start to vanish as a consequence of the hyperbolic functions present in the 
expressions. 



B. Two dimensions 



Using D = 2 in Eq. (l2~T5|) we find 

X{E) = -—j {Ecos(2nkE/u;) - sin(2irkE/uj)} . (4.10) 

■Z7T K v 71 rC ) 
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H cos(2nkE/uj) 



2tt(3 k 2 \ uj sinh 9^ 



sinh 9 k sinh 9 



cosh 9}. 



sm(2nkE/uj) } , (4.11) 



with 9% defined by Eq. ()3.34|) if we make the same approximations as in the D = 3 case. 
The ( — l) fc factor in the summand, present for D = 1, 3, is absent here, a result that is true 
for any even dimension. 

The average particle number can be shown to be 



with 

2tt CM sin(27rA;u/ti;) 7r cosh 6^ , . „ s 

N r ~ — y \- — v - 7 ; + =-1 cos(2irku/u}) \ 4.13 

fiuf^Xuj sinh0 fc /3wsinh 2 4 V P/ 7 J V ; 
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the contribution coming from the de Haas- van Alphen part of Q. 

In the continuum approximation (dropping the term in N r ) we find 

ft. ~ (2AT) 1 /2 u u \ 
to 

to leading order for large N. However the de Haas-van Alphen part of the particle number 
can be shown to lead to the step-like behaviour we saw previously in the 3-dimensional case. 
In the low temperature limit (T <C u) an asymptotic analysis of Eq. ()4.13j) can be used to 
show that 

/' 



N ~ - 
2 



ft 

La; J 



2 3 



+ 2 



-UJ. 

Solving this expression for \i gives rise to the approximate analytical solution 



UJ 



(2iV)V2 + I 



2 (4-16) 

(Again we remind the reader that the square brackets in Eqs. ()4.15|) and ()4.16|) denote the 
greatest integer function.) We therefore conclude that the step-like behaviour found earlier 
for the chemical potential when D = 3 is also found for the 2-dimensional gas. We have 
checked this result by solving Eq. (J4.12j) numerically and found that the approximation in 
Eq. (14.16)) becomes increasingly accurate as T is reduced. Because the results resemble the 
similar behaviour exhibited in Fig. ^we will not show them here. The jumps in the chemical 
potential occur, for large N, when N ~ |£ 2 for integral I. (The exact result is ^£(£ + 1) 
analogously to the 3-dimensional case studied in Ref. p.) 

It is straightforward to obtain expressions for the internal energy and specific heat. Once 
again the results are similar to those found in the 3-dimensional case, this time shown 
in Fig. |2l so need not be exhibited in detail. As T — > the result for the specific heat 
vanishes much faster than the linear approximation C ~ 7r 2 T(2A^) 1 / 2 /(3o;) predicted by the 
continuum approximation. 

V. DISCUSSION AND CONCLUSIONS 

The main conclusion of this paper is that the origin of the step-like features found in the 
numerical calculations of Schneider and Wallis P| have the same origin as the periodicity 
found in the de Haas- van Alphen effect. The general approach of Sondheimer and Wil- 

n 

son [16], so useful in the analysis of the de Haas- van Alphen effect, can be used to great 
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effect here to obtain analytical results for various thermodynamic quantities. In addition, it 
is possible to see that within this approach the continuum approximation used by Butts and 

n 

Rokhsar jj| corresponds to the neglect of an infinite set of poles in the Laplace transform of 
the partition function. It is also possible to recover the continuum approximation using the 
methods described in the general framework of Ref . 3| • 

In Sec. IIHI we obtained approximate analytical results for the 3-dimensional, isotropic 
harmonic oscillator potential. In the case T -C uj, a very simple result (Eq. (j3.21|) ) was found 
for the chemical potential that clearly exhibited the step-like features found in Ref. 5 . This 
was used to evaluate the specific heat, and it was shown (in agreement with Ref. 5|) that 
there were significant deviations from the linear temperature dependence predicted by the 
continuum approximation. The specific heat was seen to exhibit a periodic structure in the 
particle number. 

In Sec. |lV]we studied the trapped Fermi gas in 1 and 2 spatial dimensions. For the 1- 
dimensional gas we were able to show that as T —> the specific heat vanished exponentially 
fast in the inverse temperature (Eq. (J4.9)) ). The 2-dimensional gas was similar in many ways 
to the 3-dimensional case of Sec. IIHI 

Although the analysis presented above was restricted to the isotropic potential for sim- 
plicity and brevity, the same methods can be used to examine anisotropic potentials. The 
technical details are more involved due to the structure of the poles in the inverse Laplace 
transform of the partition function as mentioned above. A preliminary report of both the 
results of this paper, and the anisotropic case was given earlier j^J. In addition to a peri- 
odicity in the particle number, there can also be a periodicity when the trapping potential 
is altered. This latter effect was not found for the isotropic case and it is worth commenting 
on why this occurs, in contrast with what might be expected from the de Haas-van Alphen 
effect where the thermodynamics shows a periodicity in the inverse magnetic field strength. 
The difference between the two cases is related to the relationship between the chemical po- 
tential and the particle number. In both cases (de Haas- van Alphen and trapped Fermi gas) 
the periodic structure results from a trigonometric dependence on fi/u. (For the de Haas- 
van Alphen effect u is the cyclotron frequency associated with the magnetic field strength.) 
In the de Haas- van Alphen effect, when fi is solved for in terms of the particle number the 
leading order contribution turns out to be independent of uj. Thus the trigonometric func- 
tions that involve fx/u exhibit the familiar de Haas- van Alphen oscillations as the magnetic 
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field is varied. For the trapped Fermi gas in an isotropic potential we found /i oc oj, so that 
fi/u is independent of u, although there is a dependence on the particle number. We claim 
that this is an artifact of the simplicity of the isotropic potential, and that the situation for 
anisotropic potentials is more interesting. An extensive examination of trapped Fermi gases 
in anisotropic potentials will be given elsewhere j^]. 

The treatment presented here has only been performed for the free Fermi gas. As already 
mentioned in the introduction, the neglect of interactions for a single component gas is a 
good approximation 0, ^]. In Ref. 0] a two component model was studied that allowed for 
an interaction between the two spin components. A numerical study showed that as the 
interaction strength is increased the step-like behaviour, like that we have been discussing, 
is increasingly suppressed. For the single component gas, we expect that an analysis similar 
to that given by Luttinger 22] in the de Haas-van Alphen effect could be used to show that 
both the amplitude and period of the oscillations are not affected by interactions to leading 
order. (This conclusion does not hold for the non-oscillatory part of Q, that we have called 
fio-) It would be of interest to study the Luttinger analysis for a multi-component Fermi 
gas in more detail to make contact with the results of Ref. jf|. 
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